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1 Introduction 

The lightest supersymmetric particle (LSP) with R-parity conservation is absolutely stable and can 
contribute to the present Universe energy density as a dominant component of cold dark matter 
(CDM). If the strong CP problem is solved by introducing a light pseudoscalar, an axion, its fermionic 
SUSY partner, an axino, can be a natural candidate for CDM if it is the LSP. The relic axino CDM can 
be produced either non-thermally, through out-of-equilibrium decays, or thermally, through scatterings 
and decays in the hot plasma, as originally shown in Refs. [1, 2]. The scenario was subsequently 
extensively studied during the last decade [3-10], with either a neutralino or a stau as the next-to- 
lightest supersymmetric particle (NLSP). Ways of testing the axino CDM scenario at the LHC were 
explored in Refs. [8, 11-13] and implications for Affleck-Dine Baryogenesis in Ref. [14]. 

The strong CP problem is naturally solved by introducing a very light axion field a. The axion 
appears when the Peccei-Quinn (PQ) symmetry is broken at some scale fa- Below the PQ scale, after 
integrating out heavy quarks carrying PQ charges [15], an effective axion-gluon interaction is given 

by 

= £^«G^^G^^ (1.1) 

where as is the strong coupling constant, G is the field strength of the gluon field and G^^ = ^e^ypo-G'"' 
is its dual with £9123 = — 1- Different types of axion models have been proposed, with distinctively 
different couplings to Standard Model (SM) fields, depending on their PQ charge assignment. 
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Very light axion models contain a complex SM singlet scalar field carrying a PQ charge. In the 
Kim-Shifman-Vainstein-Zakharov (KSVZ) class of models [15] the PQ charges are assigned to new 
heavy quarks, while in the Dine-Fischler-Srednicki-Zhitnitskii (DFSZ) approach [16] the PQ charges 
are assigned to the SM quarks. This difference is the origin of different phenomenological properties [17] 
of the KSVZ and DFSZ classes of models since in the low energy effective theory at the electroweak 
(EW) scale (after integrating out heavy fields) the gluon anomaly term is the source of all interactions 
in the KSVZ models while the Yukawa couplings are the source of all interactions in the DFSZ models. 
For solving the strong CP problem, one needs a coupling of the axion to the gluon anomaly and this 
is generated by a heavy quark loop in the KVSZ models, appearing as a non-renormalizable effective 
coupling when those heavy fields are integrated out. In the DFSZ models instead the coupling is 
generated by SM quark loops. 

The couplings arising in these two popular classes of models will be discussed in the next section. 
For the KSVZ axion, constraints on the PQ scale fa have been obtained from astrophysical and 
cosmological considerations and the scale is limited to a rather narrow window 10^" GeV < /a < 
10^^ GeV [17], while for the DFSZ case precise constraints on fa have not been derived yet. 

The axino as a candidate for CDM was originally studied mostly for the SUSY version of the 
KSVZ axion model, in an important production mode corresponding to the interaction term given 
by Eq. (1.1) [2]. The supersymmetrization of axion models introduces a full axion supermultiplet A 
which contains the pseudoscalar axion a, its scalar partner saxion s, and their fermionic partner axino 
a. 



where Fa stands for an auxiliary field and i9 for a Grassmann coordinate. 

The effective axion interaction of Eq. (1.1) can easily be supersymmetrized using the superpotential 
of the axion and the vector multiplet Wa containing the gluon field 



Effective axion multiplet interactions with the other gauge bosons can be obtained in a similar way. 
Axino production from QCD scatterings due to interaction Eq. (1.3) has been considered in the 
literature in different approximations, which have led to somewhat different numerical results. The 
main technical difhculty and source of uncertainty is the question of how to regulate the infrared 
divergences due to the exchange of massless gluons. In the original study [2] a simple insertion of 
a thermal gluon mass to regulate infrared (IR) divergences was used and the leading logarithmic 
term was obtained, without much control over the subleading finite piece, since the thermal mass was 
introduced by hand. Subsequently, a hard thermal loop (HTL) resummation method was applied in 
Ref. [5] , allowing for a self-consistent determination of the gluon thermal mass and more control over 
the constant term in the high energy region of axino production. Recently, in Ref. [6] a new calculation 
was presented which, although not gauge invariant, includes more terms of the perturbative series, in 
particular the decay of the gluon whose thermal mass can be larger than the gluino and axino mass 
taken together. The two latter methods have their own advantages and limitations, and they coincide 
in the high energy region where the convergence of the perturbation series is stable. In Ref. [6] 
a previously neglected dimension-5 term in the Lagrangian, containing purely interactions between 
supersymmetric particles was also included. This however changes the axino production rate by less 
than 1%. In our paper, we adopt the original way of effective mass approximation but we improve 



—={s + ia) + V2a.'d + Fa^-&, 



(1.2) 




(1.3) 
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it to make the result positive definite. Although this method is not gauge invariant either, it is only 
known viable method at relatively low reheating temperatures, which is the regime important for the 
axino as cold dark matter candidate. We shall come back to this discussion in more detail below. 

In the calculations published so far the couplings of the axino to the gauge multiplets other than 
the one of the gluon was also neglected. In fact, in Ref. [2] a chiral transformation of the left-handed 
lepton doublets was performed to remove the axion SU{2) anomaly interaction and to leave only the 
axion U{1) contribution. Then the axion SU{2) anomaly coupling re-appears in principle from the 
leptonic loops, which are independent of the fermion masses. The corresponding axino loops, on the 
other hand, are suppressed by the ratio mf^^^^^/M'^ where M is the largest mass in the loop. Therefore 
the error in neglecting the axion SU{2) anomaly is estimated to be at most of order m^/m^ compared 
to the SU{S) term. However, in supersymmetric extensions of axionic models, the chiral rotation on 
the lepton fields involves also their scalar counterparts, the sleptons and the saxion, which generates 
additional couplings for those fields. It is also not clear if it is justified to perform a redefinition in 
a fully supersymmetric way when supersymmetry is in any case broken. Rather than going into the 
swamp of such interactions with unknown slepton and saxion masses as parameters, we will keep here 
the general axion SU{2) anomaly interaction and its SUSY counterpart, which is more tractable, and 
examine its effect on the axino abundance. 

Moreover, in the present work we will also consider the role of Yukawa-type interactions between 
the axino and the matter multiplets, which arise either at one loop or at tree level, depending on the 
model. We will investigate what effect these terms may have and in particular how model dependent 
our results for the axino abundance are. In the previous study, only the KSVZ model was considered 
for the axino production. However in the DFSZ model axino production is dominated by the Yukawa 
coupling and the dependence on the reheating temperature is quite different from that in the KSVZ 
model. Only recently Refs. [33] and [47] considered axinos in the DFSZ model. In particular, Ref. [33] 
gave a simple approximate formulae for the relic density of light axino as dark matter. In our paper, 
we study axino production in the DFSZ model and the suitability of the DFSZ axino as dark matter 
in a more complete way. 

In view of the recent developments in estimating the QCD contribution, in this paper we also 
update and re-examine relic CDM axino production, with an emphasis on an estimate of the un- 
certainties as well as on model dependence. In particular, our updated analysis of the axino CDM 
scenario improves on the previous works in the following aspects: 

(i) an inclusion of some previously neglected terms in the axino production and of subleading terms 
in the mass of the gluon beyond the logarithmic and constant pieces - these last parts ensure that the 
cross-section remains positive even for the invariant mass s smaller than the gluon thermal mass; 

(ii) an explicit calculation of axino production via SU{2) and U{\) interactions; 

(iii) a derivation of the axino abundance in specific implementations of the KSVZ and DFSZ models. 

(iv) an update on the constraints on the reheating temperature Tp. ^ for the both the neutralino and 
the stau as the NLSP using the current WMAP-7 result on the DM relic density and relevant structure 
formation data. 

In Sect. 2, we define the axino by its relation to the shift symmetry of the KSVZ and DFSZ axions. 
In Sect. 3 we calculate axino production rate for the KSVZ and DFSZ axion models, and in Sect. 4 

^ We assume the instant reheating approximation and define the reheating temperature as the maximum temperature 
at which standard Big Bang expansion with a thermalised bath of SM particles starts. We can easily translate the axino 
abundance given with this conventional definition of the reheating temperature to more specific reheating scenarios. 
E.g., Ref. [6] considers a reheating process from the decay of the inflaton and obtains a slightly smaller abundance of 
axinos, reduced by a factor of 0.745. In comparing our results we account for this difference. 
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we discuss several constraints on the scenario arising from cosmology. Finally, in Sect. 5 we present 
our conclusions. 



2 Axinos 



2.1 The framework 



In a recent review [17] low energy axion interactions were given in terms of the effective couplings 
with the SM fields ci, C2, and C3, which arise after integrating out all heavy PQ charged fields. The 
effective axion Lagrangian terms are 



fa 

T~ X] {c2'^u^i'i75Ui + d^m'-^dii-f^di) (2.1) 



fa 

C3 



327r2/a 



aGG, 



where C3 can be defined to be I (if it is non-zero) by rescaling fa-^ The ci term is the PQ symmetry 
preserving derivative interaction of the axion field that can be reabsorbed into the C2 term by a 
partial integration over on-shell quarks. For the C2 terms, we have only kept the lowest order terms 
(proportional to 1/ fa), while in principle an infinite series of terms in a/ fa arises. 

In the following we will consider the two popular scenarios mentioned earlier: the KSVZ and 
the DFSZ classes of axion models. The KSVZ class of axion models corresponds to the choice ci — 
C2 = 0,C3 = 1, and the DFSZ one to ci = C3 = and C2 ^ 0, after integrating out the heavy field 
sector responsible for PQ symmetry breaking. In the latter model, if the Higgs doublets Hu^d carry 
respective PQ charges Qu,d, the SM fields also carry PQ charges,'^ see Table I, and the anomaly 
interaction proportional to C3 ^ arises from SM quark loops. The axion mass is given by the strong 
interaction and is proportional to |c2 + C3I. Hence the sum C2 + C3, if it is non-zero, defines the QCD 
axion. Only this combination of the two couplings is physical, since a chiral axion-dependent PQ 
rotation of the quark fields can shift the values of C2 and C3, while keeping C2 -I- C3 constant. This is 
connected to the well-known fact that, if one of the quark masses is zero then the anomaly becomes 
unphysical and can be reabsorbed in the rotation of the massless field. 

In the supersymmetric version of axionic models, the interactions of the saxion and the axino with 
matter are related by supersymmetry to those of the axion. Hence the definition of the axion at low 
energy must be connected to the definition of the axion multiplet, and therefore of the saxion and the 
axino, at energies above the EW scale.* 

Below we will examine more closely the KSVZ and DFSZ models of the axion. In both models 
one imposes the PQ symmetry at a high energy scale and the axion emerges from its spontaneous 
symmetry breaking. As a specific example, let us consider the PQ sector at high energy with the PQ 



^The strong coupling 33 is usually omitted by absorbing it to the field strengths. Here, fa denotes the so-ealled axion 
decay constant which is typically related to the vacuum expectation value V of the PQ symmetry breaking scalar field 
by fa = V/N where N is the domain wall number. 

^In variant axion models, Q„ ^ may have family dependence but here we suppress family indices. They can be 
inserted if needed. 

*The supersymmetrization of axion models was first discussed in Refs. [18-20]. An explicit model was first constructed 
in [21], and the first cosmological study was performed in Ref. [22]. 
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Model 


Z 


Si 


S2 


Ql 


Qr 




Hu 




^R 


U'r 


KSVZ 





Qa 


-Qa 





















DFSZ 





Qa 










-Qd 


Qu 





+Qd 


+Qu 



Table 1. The PQ charge assignment Q. Ql and Qj^ denote new heavy quark multiplets. In the DFSZ model 
2Qct = Qu + Qd, and the PQ charge of the left-handed SM quark doublets ql vanishes. See text for more 
details. 



breaking implemented as in Ref. [21] 



WpQ = fzZ{SvS2 - K'), 



(2.2) 



where Z, S'l, and S2 are gauge singlet chiral superfields, fz is a Yukawa coupling and Va is a parameter 
in the Lagrangian which determines non-zero VEVs (S'l) = (6*2) = Va in the minimization process. 
The superfields transform under the J7(1)pq symmetry as 



Si 



'Si, 5*2 



(2.3) 



The potential in the global supersymmetric limit is 



dW 



d(pa 



1 



-D^D" 



(2.4) 



with 



(2.5) 



where g is the gauge coupling of the gauge groups and 
the superpotential given by Eq. (2.2) one has 



' denotes collectively all the scalar fields. With 



dW 
'dZ 
dW 

dS^ 

-ds;-^'^^'- 



= fzZS2, 



(2.6) 



At the tree-level both S'l and S2 develop VEVs and break the PQ symmetry. The fermionic partners 
of Z and S' = [Si + S2)/\/2 combine to become a Dirac fermion with mass mz = V^fzVa, while 
S — {Si — S2)/V2 contains the axion field and can be identified with the axion multiplet. From 
Eq. (2.6) we note that, if {Z) = and SUSY is not broken then both the axino and the saxion are 
mass degenerate with the axion. 

However, with soft SUSY breaking terms included, (Z) can develop a non-zero value, in which 
case both the saxion and the axino become massive, independently of the axion. Therefore, a full 
specification of the SUSY breaking mechanism is needed in order to determine the mass of the axino 
exactly. This was first studied in Ref. [21] for the superpotential above, while another superpotential 
and SUSY breaking with a very small axino mass was given in Ref. [23] . Recently the case of a direct 
coupling between the axion and SUSY breaking sector was discussed in Ref. [24]. 
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2.2 The KSVZ model 

In the KSVZ approach in order to obtain the anomalous interaction of the axion and the gluon fields, 
one introduces the heavy quark fields Ql and Qe. in the superpotential as Ref. [15] 

M^KSvz = W^PQ + IqQlQrSi. (2.7) 

After the [/(1)pq symmetry breaking takes place, as discussed above, Ql and Qr combine to become 
a heavy Dirac fermion with mass mq — fqVa, Assuming that Va 'S> ^13/2, the scalar partner has 
practically the same mass and the whole supermultiplet can be integrated out in a supersymmetric 
way. The low-energy Lagrangian can then be obtained by integrating out the heavy quark multiplet 
or by using anomaly matching condition. In this way one finds that the axion anomaly term becomes 

r-ff = ^£^aG^''G^., (2.8) 

where the axion a is the pseudoscalar component of the superfield S = {Si — S'2)/-\/2, and the axion 
decay constant is fa = 2t4. Nq is the number of the heavy quarks and we consider Nq = 1 in our 
case. Then we have ci = C2 = and C3 = 1 [17]. 

The low-energy interactions of the saxion and the axino fields can be obtained in the same way 
by integrating out the heavy (s)quark fields. However, in the limit of unbroken SUSY the low-energy 
effective Lagrangian should be given in a SUSY invariant form, including Eq. (2.8), as 

^off ^ ( ^2 p^aa^a ^ ^ ^2.9) 



The axion superfield A and VF" " are given by [25] 
The effective Lagrangian in terms of the Bjorkcn-Drell gamma matrices then reads 



(2.10) 



cff 



a(G'^^^G^, + i5,,(g7V.9)) 



I i7.r«^ [t'^jT ]_,5~ 
+^aG^^ 7 g- 

+ \/2(FAAA + h.c.)., 
where the giuino and axino 4-spinors are given by 



(2.11) 



9=[ f), (2.12) 



The effective Lagrangian (2.11), including the axino interaction with the D-term has recently been 
used in Ref. [6] to calculate the thermal production rate of axino dark matter and will be the basis 
also for the present analysis. 



^ Our saxion-gluino-gluino interaction differs by a factor of 2 from that in Ref. [6] . 
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2.3 The DFSZ model 



In the DFSZ framework, the SU{2)l x U{1)y Higgs doublets carry PQ charges and thus the Ught 
quarks are also charged under C/(1)pq. The charge assignment is shown in Table 1. The anomaly 
coupling aGG can be obtained after electroweak symmetry breaking (EWSB) through the coupling 
of the axion to the Higgs doublets which couple to the light quarks. To this end, one adds a non- 
renormalizable term to the PQ breaking superpotential [26] 



W^DFSZ = W^PQ 



A. 

Ah 



■SiHdHu, 



where Wpq is given in Eq. (2.2) and HdHu = eapH^H!^. Note that here fj.Mp/V^ 
a phenomenologically acceptable supersymmetric yit-term. 
The superpotential including light quarks is given by 

VKmssm = ytQHuU' + ybQHdD' + VrLHaE'. 



(2.13) 

1 generates 

(2.14) 



Before the EW symmetry is broken but after the J7(l)pQ symmetry is broken, the massless axion 
Goldstone boson is identified as 



a — -= — . 

Here we defined the component fields in the same way as in Eq. (2.10), 



(2.15) 



Si 



V2 



(2.16) 



and similarly for Instead, the axino mass eigenstate can be obtained from the mass matrix 





' 




(^2)] 




"0 


zq 


Va 


fz 









-fz 


zq 





Va 




[{S2) 


(^i> 







Va 


Va 


0. 



(2.17) 



in the basis of (ai, 02, Z) and we have used zq — {Z). The lightest state, which we will identify with 
the axino, has mass fzzo and is given by 



ai — 02 



V2 



(2.18) 



which coincides with Eq. (2.15). Since EW symmetry is not broken, the axion (and thus the axino) 
do not mix with the Higgs (and higgsino) and therefore the axion cannot have the SU{?i)c anomalous 
interaction with gluons generated at one loop via the quark triangle diagrams.^ However, SU{2)l 
and U{1)y anomalous interaction can be generated via a higgsino triangle loop through the Yukawa 
coupling derived from Eq. (2.13). They are given by 



C 



anomaly 



8^ 



(2.19) 



These anomaly interactions appear also for the axino, but may obtain corrections of order of one from 
SUSY breaking, compared to the axion couplings, and such uncertainties will be later encoded in the 
coefficients Gaww ^GaVY ■ 



®An axino-gluino-gluon coupling can arise at two loops through the Higgs Yukawa couplings, which are non-vanishing 
even above EW symmetry breaking, but it is strongly suppressed and we will neglect it here. 
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After the EW and the J7(1)pq symmetries are both broken, the axion is identified with the 
Goldstone boson of the broken C/(1)pq symmetry, given by 



_ 2QaVaas - QdVdPd - QuVuPu 

where Ug = (ai — a2)/\/2, and we expanded the Higgs field as 



(2.20) 



h: 



_ Vd + hd iP^I^ 

■ V2 



V2 



(2.21) 



with Vd/^/2 ~ (0|i?(i|0) and u„/a/2 = (0|77„|0), while Pd,u are the pseudoscalar fields contained in the 
electrically neutral component of Hd,u- 

The neutral Higgs boson component eaten by the Z-boson and the orthogonal pseudoscalar Higgs, 
i.e. the phases of Hd, and Si — 5*2, are given by 

-VdPd + VuPu 



A = 



VdVu^s + VaVuPd + faVdPu 



_ + ^)Vaas - VuPd - VdPu 

_ v'^Vgas - VdvlPd - VuVdPu 

Equating (2.20) with the last term of Eq. (2.22), we obtain 



(2.22) 



Qd= — = tan/3, Q 

Vd 

2^Vu^ Vd^ 2sin/3cos/3' 



Vd_ 
1 



1/tan^, 



(2.23) 



up to a common normalization constant. 

The interactions of the axion with the matter fields are obtained through the axion part of the 
phase of the Higgs fields. 



Pd 



Pu 



= — sin /3 cos (3 + ■ 

V o 

= — sin /? cos l3 + 



Considering the Yukawa interaction from the superpotential, Eq. (2.14), 

C = -ytURULHl - yluLUuHl* H , 

the Lagrangian terms for the up-type quark axion couplings are given by 



C 



yt-^URUL exp 



V^Va 

_.w|_a_ 



(2.24) 



(2.25) 



(2.26) 



+ 
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and similarly for the down-type quarks. We can now compare the above to the coefficients in the 
definition of the axion [17] and obtain for following values for €2''^; compare Eq. (2.1), 



= ^ = cos^ /3, 4^^ = sin^ (3. (2.27) 

V w 



After integrating out all the Higgs fields except the axion supermultiplet, which remains light, all 
the axion couplings arise from the C2 terms given in Eq. (2.26) and Eq. (2.27) at low energies above 
the quark masses. At one loop in the SM fcrmions one obtains the axion-anomaly interaction term. 
It is then given by 

■^anomaly — 7^ 77vr7~T^^^ ^ fiu i (2.28) 
»7r(2\/a) 

where iV = 6 and again fa — 2Va/N . 

In the supersymmetric limit, below EW symmetry breaking, the axino mass eigenstatc can be 
read off from Eq. (2.20) and is given by 

Here h^^u denote the fermionic components of the electrically neutral parts of H^^u- However since 
supersymmetry is broken, in general in the DFSZ models the axino mixes with the higgsinos differently 
than the axion with the Higgs and the mass eigenstate is not exactly the state given in Eq. (2.29). 
Such a field is a good approximation to the physical axino only if the mixing generated from the 
superpotcntial (2.13) and SUSY breaking, which is of order Vu,d/ fa and z^/ fa, is negligible. Otherwise 
the whole axino-neutralino mixing matrix has to be considered in detail. We will not discuss this case 
further, but point instead to the related studies in the case of the Next-to-Minimal Supersymmetric 
Standard Model with a singlino LSP [27]. 

In the DFSZ models, there are also axino tree-level Yukawa interaction terms to the quark and 
squark with a coupling of the order of mg/ fa below the EWSB scale, with the Higgs and higgsino with 
coupling /i//a even above the EWSB scale. These tree- level interactions are not suppressed by a gauge 
coupling or a loop factor, as the QCD anomaly term has, and thus they give the dominant contribution 
to axino production through the decay and/or scattering processes at low reheating temperature Tr. 
At high reheating temperatures above EWSB instead the SU{3)c anomaly coupling is absent and the 
SU(2)]^ anomalous interaction dominates the production. 



3 The production of axinos 

There are two efficient and robust ways of populating the early Universe with axinos. Firstly, they 
can be produced through scatterings and decays of particles in thermal equilibrium. This mechanism, 
which we call thermal production (TP) depends on the reheating temperature after inffation. The 
other mechanism, which is independent of Tr, involves non-thermal production (NTP) of axinos from 
the decay of the NLSP after it has frozen out from the plasma. Note also that, even though squarks 
are normally not the NLSP and remain in thermal equilibrium, for Tr < ruq and large gluino mass, 
axino yield from decay processes q qa can dominate the abundance [3] . Additionally the decay of 
inflaton or moduli can produce axinos but such contributions are very model dependent and won't be 
considered here. 
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Thermally produced axinos in the keV mass range were considered as warm dark matter (WDM) 
in Ref. [28] and much lighter ones as hot DM in Ref. [22]. In Ref. [1] it was shown that axinos from 
neutralino NLSP decays can be a natural candidate for CDM and this was extended in Ref. [2] to 
include TP. If the axino mass is between around an MeV to several GeV, the correct axino CDM 
density is obtained when is less than about 5 x 10^ GeV [2]. At higher Tr and lower keV) mass, 
axinos could constitute WDM. As a digression, we note that, when axinos are very heavy, and it is the 
neutralinos that play the role of the LSP, their population from heavy axino decays could constitute 
CDM [29], leading in particular to the possibility of TeV-scale cosmic ray positrons, as pointed out in 
Ref. [30]. The possibility of either WDM or very heavy axinos is not discussed in this paper. 

The interactions leading to CDM axinos were extensively studied in terms of cosmological impli- 
cations in Refs. [2-4] and collider signatures in Refs. [8, 11-13]. If the LHC does not confirm the decay 
of heavy squarks or gluino to a lighter neutralino, the axino CDM idea with the R-parity conservation 
can not be saved unless some other mechanisms such as an effective SUSY is introduced [31]. 

In general the couplings of the axino to gauge and matter fields are analogous to those give in 



Here we normalize the PQ scale by 2Va/A^ fa which sets — 1 for the QCD anomaly coupling and 
therefore defines the coefficients Caww and CaVY in the axion models considered above. 

3.1 Thermal production 

At sufficiently high temperatures (> 10^ GeV) axinos can reach thermal equilibrium with SM particles 
and their superpartners. However, assuming that a subsequent period of cosmic infiation dilutes the 
population of such primordial axinos (and that they are not produced directly in inflaton decay), a post- 
inflationary axino population comes firstly from the hot thermal bath. If the reheating temperature is 
very high, above the decoupling temperature of axinos (~ 10^ GeV), their relic number density reaches 
again thermal equilibrium and is the same as that of photons. In that case axinos must be so light ( 
< 1 keV) that they become warm or hot DM [28]. On the other hand, when the reheating temperature 
is below the decoupling temperature, axino number density is much smaller than that of photons and 
its time evolution is well described by the Boltzmann equation without backreaction [2]. 

In the KSVZ model, at high temperatures, the most important contributions come from two-body 
scatterings of colored particles into an axino and other particles. At lower reheating temperatures, on 
the other hand, the decay of squarks or gluinos can dominate the production of axinos [3, 4]. In Ref. [2], 
an effective gluon thermal mass was introduced to regulate the infrared divergence in the scattering 
cross-section. Subsequently, a more consistent calculation using the hard thermal loop (HTL) approx- 
imation was applied to the axino production in Ref. [5]. However, the HTL approximation is valid 
only for small gauge coupling, g <C 1, which corresponds to the reheating temperature Tr >• 10^ GeV. 
Below 10^ GeV the HTL approximation becomes less reliable [5]. In fact, the production rate becomes 
even negative at gs > 1.2, and therefore the result becomes unphysical. Strumia tried to improve this 
behaviour by using the full resummed finite-temperature propagators for gluons and gluinos in the 



Eq. (2.9) 




(3.1) 
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loop [6]. This procedure includes axino production via gluon decays, which is kinematically allowed 
by thermal mass, and results in an enhancement compared to the HTL approximation. However, his 
method is gauge-dependent in the next-to- leading order [32], indicating that not all the contributions 
of that order are included, and therefore also does not give a fully satisfactory result in the large 
gauge coupling regime. For these reasons, we believe that it is worth pursuing an alternative method 
of computing the rate at large couplings and we apply the effective mass approximation as the only 
way to evaluate the axino thermal production at large coupling g and low reheating temperature 
^ 10^ GeV. Even though this method is also not gauge invariant, it captures relevant physical 
effects like plasma screening and gives positive and physical cross-sections for each single scattering 
process. 

As stated earlier, in the previous studies of TP of relic axinos the contributions from SU{2)l 
and U{1)y gauge interactions were neglected, but for completeness, we will discuss here all SM gauge 
groups explicitly in order to examine their possible effects. 

To start with, in evaluating the contributions from the strong interactions we will follow the 
method used previously in Ref. [2] to obtain the (dominant) axino TP cross-section. We will further 
update and correct the tables presented there by following Ref. [6] to include the previously ignored 
dimension-5 axino-gluino-squark-squark interaction term. This term changes only the contribution 
from the processes H and J in Table 2, while it does not affect the other terms. The processes B, F, 
G and H, with gluon t- or w-channel exchange, are infrared- divergent and thus, following Ref. [2], are 
regularized with the inclusion of an effective gluon thermal mass in the gluon propagator. This method 
gives always positive definite values for the single cross sections. The full results and the method how 
they are obtained is described in the Appendix. In Table 2, we give for simplicity only the first two 
leading terms in the expansion for s/m^g >> 1. However the logarithm in the approximate formulae 
of Table 2 gives unphysical negative value for s < rri^g. Therefore in the numerical calculation we 
have used the full formulae which are positive definite for all values of s listed in the Appendix. 

The total cross-sections cr„, where n = A, . . . , J , can be written as 

where 5-„(s) denotes the cross-section averaged over spins in the initial state and are given in Table 2. 
The relevant multiplication factors are also listed: rtg (the number of initial spin states), np (the 
number of chiral multiplets) and rji (the number density factor, which is 1 for bosons and 3/4 for 
fermions). We assume particles in thermal equilibrium to have a (nearly) Maxwell-Boltzmann distri- 
bution, proportional to r]i, and neglect Fermi blocking or Bose-Einstein enhancement factors, which 
are close to one at these temperatures. We restrict ourselves to temperatures above the freeze-out 
temperature of SM superpartners involved, so that the approximation of thermal equilibrium is always 
satisfied.^ Finally, in Table 2 the group theory factors /"'"^ and Tjjj, of the gauge group SU{N) satisfy 

the relations Ea.f,,c If'^l' = ^(^' " 1) and Ejk l^/J' = " l)/2- 

Next we move to include contributions from the SU{2)l and U{1)y gauge interactions. The 
relevant axino-gaugino-gauge boson and axino-gaugino-sfermion-sfermion interaction terms, in view 

^At lower temperatures, superpartner number densities, apart from the NLSP one, drop down to zero and the NTP 
regime is reaehed. 
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n 


Process 








V1V2 


A 


+ ^ a + r 




4 


1 


1 


B 


g°- ^ a + g" 




4 


1 


3 
4 


C 


g"^ + qk^a + q.j 




2 


2Np 


1 


D 


g°- + qk^a + qj 


1 \rpa |2 
32 l-'ifcl 


4 


Nf 


3 
4 


E 


qj +qk^a + g" 


J^\rra 12 


2 


Np 


3 
4 


F 


5" +g'' ^ d + g" 




4 


1 


3 3 

4 4 


G 


g'' + qk ^ a + qj 


i|T;^p[log(s/m2^)-2] 


4 


TVi. 


3 3 

4 4 


H 


9°^ + qk^a + q.j 


i|T;5^|2[log(./m2J-|]* 


2 


2Nf 


3 
4 


I 


qk + qj ^ a + .9° 




4 


^Nf 


3 3 

4 4 


J 


qk + q* -^5, + g" 




1 


Nf 


1 



Table 2. The cross-section for each axino thermal-production channel involving strong interactions. The 
particle masses are neglected except for the plasmon mass moff. The H and J entries with an asterisk in the 
third column are changed due to including the missing term and cross-sections or n_F in the others processes 
(A,B,D,E,F, and I) are corrected from those of Ref. [2]. The logarithm in the approximate formulae in this 
Table gives unphysical negative value for s < nr^g. Therefore in the numerical calculation we used the full 
formulae which are positive definite for all values of s. The full cross sections for the processes B, F, G and H 
are given in the Appendix. 



of Eq. (2.11), are given by 



167r/a ^ 



-^^^'E^JhT^fo (3.3) 

fn 

+ » 1^ f 0757^,7 



where the terms proportional to a2 come from the SU (2) f and the ones proportional to ay from the 
U{1)y gauge groups. Caww and CaYY are the model-dependent couplings for the SU{2)f and U{1)y 
gauge group axino-gaugino-gauge boson anomaly interactions, which is defined after the standard 
normalization of /a, as in the first line for S'[/(3), as stated below Eq. (3.1). Here 0:2, W^u and 
ay, F, Y^i, are the gauge coupling, the gaugino field and the field strength of SU{2)l and U{1)y 
gauge groups, respectively, fo represents the sfermions of the iS'J7(2)-doublet and / are the sfermions 
carrying the U{1)y charge. 

We start from Table 2 and replace quark triplets with S'?7(2)L-doublcts with corresponding group 
factors. For the abelian U{1)y factor, the processes A, B, and F vanish and we can replace \T°-f.\'^ 
with the square of the corresponding hypercharges. Finally, we use Np = (12,14,11) to count the 
matter multiplets charged under the MSSM gauge groups {SU{3), SU{2)f and U{1)y, respectively. 
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We include above the SUSY breaking scale the full 1-loop MSSM running of the gauge couplings and 
gaugino masses. 

The second term for each gauge group in Eq. (3.3) also generates three-body gaugino decays into 
an axino and two sfermions, assuming that the gauginos are heavy enough. The three-body decay rate 
of the gluino is given by 



r(g° ^ aq,q*k) 



where nig denotes the gluino mass, 



1^ 128^3 J2 j[iQT,\ 



(3.4) 



G{x) =-\/l -4a;(l + 5x - 6x^) 
3 

- 4a;(l - 2.T + 2x^) log 



1 -f yi - 4x 
2x 



- 2 



(3.5) 



and the mass of the axino has been neglected. However, the three-body decay is suppressed by an 
additional power of the gauge coupling constant and is kinematically allowed only when the gluino 
mass is larger than the sum of the two final-state squark masses. Therefore gluino three-body decay 
through the second term in Eq. (3.3) is subdominant to the two-body decay. 

As stated in Ref . [3] , an effective dimension-4 axino-quark-squark coupling can be generated at a 
loop level also in the KSVZ model. Here we take into account this effective Yukawa interaction, which 
appears at a two-loop level in the KSVZ models and a tree level (with a tiny mixing) in the DFSZ 
models [3], 



(3.6) 



where ijjj and i/^j denote the SM fermions and their superpartners. 

In the KSVZ class of models, the effective coupling comes predominantly from the logarithmically 
divergent part of the gluon-gluino-quark loop term and is proportional to mg [3] , 



L/R 
5cff 



V27r2 fa 



9 1 I fa 

"log — 



(3.7) 



Subleading terms have not yet been computed and may give a correction of the order of 20 — 30%, in 
analogy with what has recently been obtained in Ref. [8] for the effective tau-stau-axino coupling. 

In the DFSZ models there exists also a tree-level axino-quark-squark coupling which is proportional 
to the mass of the quark [33], coming from the C2 interaction term in Eq. (2.27), 



L/R. 
9eS 



fa 



cos^ f3 
fa isin^/?, 



(3.8) 



where the upper row relates to the up-type quarks and the lower row to the down-type quarks. These 
tree-level couplings are always smaller than the one-loop ones for the light generations, but not for 
the third one. In fact for the top quark, the tree-level coupling dominates if tan /3 < 4 for the gluino 
mass of 700 GeV, while the bottom quark tree-level coupling dominates if tan /3 > 1 for the same 
choice of the gluino mass. Note that, at low reheating temperatures, only the top-stop-axino coupling 
is important for axino thermal production. This is because the lighter stop is usually the lightest 
colored superpartner and remains in equilibrium to rather low temperatures, even below the EWSB 
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scale. Similarly, there exists an effective tau-stau-axino vertex, which was first obtained in Ref. [4] 
and more recently re-derived in Ref. [8] in a full two-loop computation. This coupling is smaller and 
not important for thermal production, but it is important for the non-thermal production when the 
stau is the NLSP. 

In the DFSZ models, there is also a tree-level axino-Higgs-higgsino coupling [47]; compare the 
second term in Eq. (2.13), 

« = (") 

where fi ^ fsV^/Mp and fa = 214- It contributes to axino TP through higgsino decays in thermal 
equilibrium, and can be comparable to that from squark decays through Eq. (3.7) and Eq. (3.8), or 
even become much larger if /i is larger than the top quark mass. The axino production due to this 
coupling in DFSZ models has been recently investigated also in Ref. [33]. Note that the Yukawa 
coupling and the axino-higgsino mixing contained in the neutralino mass matrix also give rise to 
additional scattering channels contributing to the axino production, but we will neglect them here 
since such dimension-4 scatterings are usually less important than the decays [3]. 

We have evaluated the thermal production of axinos numerically and present the results in Fig- 
ures 1 and 2 for representative values of fa = 10^^ GeV and rriq — rrig = 1 TeV. We do not consider 
here the dependence on masses, however see Ref. [3]. For different values of fa the curves move up 
or down proportional to In Fig. 1, we show the axino yield Y (where Y = n/s is the ratio of 
the number density to the entropy density) from strong interaction in the KSVZ model. Our result 
obtained with the effective mass approximation is shown with the solid black line. Compared to the 
previous plot in Ref. [2], the inclusion of the squark decay changes the plot at low reheating temper- 
ature, while the other new squark interactions do not have any noticeable effect. There is a factor 3 
difference in the abundance at high reheating temperature compared to that in Figure 2 of Ref. [2], 
which was a numerical error at that time and was corrected later. For comparison, the axino yield 
from scatterings using the HTL approximation [5] is plotted with the blue (dashed) line and Strumia's 
result [6] is shown with the green line. 

For Tr > 10'' GeV the axino abundance using the effective mass approximation increases consis- 
tently with that of Strumia. We found that the difference between the two prescriptions is of order a 
factor of three. In principle we could reabsorb this difference in the definition of the effective gluon 
mass at high temperature, which in our scheme cannot be determined self-consistently. With this 
tuning we could match the perturbativc result at high temperature. On the other hand, doing this 
would require a gluon thermal mass smaller than the expression ~ gT used in our calculations. Hence 
we prefer instead to consider this factor as an estimate of the theoretical error of using the effective 
mass approximation at high reheating temperatures. We assume that this error does not increase for 
reheating temperatures less than 10* GeV and we apply the effective mass approximation to the DFSZ 
model. 

For lower temperatures the contributions from the decays of squarks and gluinos in thermal plasma, 
which were not included in Ref. [6] , start playing some role. We mark those in Fig. 1 with a green solid 
and red dashed curves, respectively. It is known that, at reheating temperatures above superpartner 
masses scattering diagrams involving dimension-4 operators are usually subdominant relative to those 
coming from dimension-5 operators and to decay terms, and are negligible. Also the decays do not 
give significant contribution to the TP of axinos, apart from very low Tr [3], and this is confirmed in 
Fig. 1. While all the above contributions are generated by strong interactions only, for comparison, 
we show also (as a black dashed line) the relative contribution from an out-of-equilibrium bino-like 
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Figure 1. Thermal axino yield Yg as a function of the reheating temperature Tr from strong interactions 
using the effective mass approximation (black). We use the representative values of fa = lO'^^ GeV and 
rriq — rrig = ITeV. For comparison, we also show the HTL approximation (dotted blue/dark grey) and 
that of Strumia (green/light grey). We also denote the yield from squark (solid green/light grey) and gluino 
decay (dotted red), as well as out-of-equilibrium bino-like neutralino decay (dashed black). Here we used 
the interactions in Eq. (3.3) and Eq. (3.7) for the KSVZ model. We use the same definition of reheating 
temperature in the instantaneous reheating approximation for the three methods. 



neutralino decay to an axino and a photon originally considered in Ref. [1]. It is clear that NTP is 
only important at very low Tr, well below squark or gluino masses. 

In Fig. 2, we show a contribution to the axino yield in thermal production from each SM gauge 
group interaction. Here we set the coefficients Caww = 1 and CaVY = 1 as a normalization. As 
shown in the figure, the contributions from scatterings due to SU{2)l and U{1)y couplings (blue 
dotted and green solid lines, respectively) are significantly suppressed compared to that of SU{3)c 
(red dashed), by a factor of 10 or more. This is because the interaction between axinos and gauge 
bosons are proportional to a gauge coupling-squared so that the cross-section is cr cx a^. Thus it 
would be only for very large (and perhaps unnatural) values of the effective couplings C'aww j CaVY 
that these channels could become comparable to the QCD contribution. To give an order of magnitude 
estimate of these effects, we included the SU(2) and U(l) contributions with a normalized value in 
figure 2. For different values of CaWW and CaYY the curves move up and down. We note that in 
general SUSY breaking effects in the leptonic sector may bring a modification of the couplings here 
considered. The situation here is different from the case of gravitino production since the interactions 
of the gravitino to the three gauge groups are of the same order; the spin- 3/ 2 gravitino component 
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Figure 2. Thermal axino yield Yg as a function of Tr, from each of the SM gauge groups. Here, we have 
used C'aww = C'aYY = 1- The lines at high Tr are not perfectly parallel due to the running of the gauge 
couplings, which affects the SU(3) yield more strongly and in the opposite direction than the other gauge 
groups. 



couples in fact universally, while the goldstino component proportionally to the gaugino masses. We 
therefore conclude that the QCD contribution is strongly dominant in the KVSZ models and so the 
axino production at high Tr is practically model independent as long as the number of heavy PQ 
charged states can be absorbed into the definition of fa- 

However, at low Tr the thermal production from scatterings becomes strongly suppressed by the 
Boltzmann factor. In the region where Tr < 100 — 1000 GeV, axino production due to the decays of 
gaugino, squarks or neutralinos become important. Actually, the lightest neutralino decay via U{1)y 
couplings becomes dominant in the very low reheating temperature regime since the number density 
of the heavier colored particles becomes strongly suppressed by the Boltzmann factor there. On the 
other hand the neutralinos, depending on their composition and the supcrsymmetric spectrum, can 
freeze-out with a still substantial number and then give rise also to non-thermal axino production, as 
we have seen in Fig. 1. This contribution to the axino yield is usually more important than the one 
due to neutralino decays in equilibrium, which is proportional to CaVY and typically below lO"'^^. 

For the case of the DFSZ instead, the role of the QCD interaction is played by the SU{2) interaction 
and the dominant decay term above the EW symmetry breaking is the higgsino one instead of the 
squark one. 

Our results for the total thermal production yield for both KVSZ and DFSZ type of models can be 
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Figure 3. Thermal axino yield Y^^ as a function of Tr for two specific KSVZ models: Qcm ~ (CuYY = 0) 
and Qem = 2/3 {CaYY = 8/3), and for a DFSZ model with the {d'',e) unification [34], for which we used 
/i — 200 GeV and the higgsino mass mf^ — 200 GeV. The horizontal lines show the values of axino mass for 
which the corresponding axino abundance gives the correct DM relic density. 

seen in Figure 3.^ There we show the KSVZ model with different values of the CaYY couphng. In the 
case of non-zero CaYY the contribution from neutrahno decay in equilibrium can be clearly seen for 
Tfl ^ 10 GeV and it is very suppressed. Moreover we give also the yield for the DFSZ model in solid 
green, for /i = 200 GeV and the higgsino mass — 200 GeV. We can see that even for relatively small 
fj,, axino production from higgsino decay dominates over the one from the anomaly terms for reheating 
temperature Tr < 10^ GeV. The importance of axino-Higgsino-Higgss coupling in the DFSZ model 
was recently discussed in [33] and our result is consistent with that analysis. The abundance is so 
large that the CDM density can be reached with an axino mass as small as 100 keV, independently of 
the reheating temperature.^ In this range of reheating temperature, the axino production from decay 
dominates that from scatterings. Therefore the use of the effective mass approximation or another 
IR screening prescription in the scattering process is irrelevant to the axino production in the DFSZ 
model in the range of reheating temperature where the decay dominates. For higher Tr, the SU{2)l 
anomaly term starts dominating and the abundance is proportional to Tr as in the KVSZ case, but 
with a smaller coefficient. In the same figure, we also mark horizontal lines corresponding to the axino 

®A similar figure for the DFSZ model is given in Ref. [47]. 

® Such general effect due to decaying particles in equilibrium has been recently called "freeze-in" in Ref. [35] and 
discussed for the axino in Ref. [36] . The freeze-in mechanism was included already in gluino or squark decays to axinos 
in the plasma in Ref. [3] . 
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Figure 4. Tr versus ma for fa = 10" GeV in the KSVZ models. The bands inside like curves correspond 
to a correct relic density of DM axino with both TP and NTP included. To parametrize the non-thermal 
production of axinos we used Inlsp = (1), 10"^" (II), and 10~* (HI)- The upper right-hand area of the plot 
is excluded because of the overabundance of axinos. The regions disallowed by structure formation are marked 
with vertical blue dashed lines and arrows for, respectively, TP (ma < 5keV, see text below Eq. (4.2)) and 
NTP (ma < 30MeV, with a neutralino NLSP). 

mass giving the correct DM relic density for the given relic abundance of Fa • 
3.2 Non-thermal production 

As stated above, axinos can be produced non-thermally in NLSP decays after they have frozen out 
of equilibrium. This NTP mechanism is dominant for reheating temperatures below the mass of the 
giuino and squarks [1,2]. In this case, the axino abundance is independent of the reheating temperature 
as long as the temperature is high enough for the NLSP to thermalize before freeze-out. Axino relic 
density from NTP is simply given by 



Clearly, in order to produce a substantial NTP population of axinos, the NLSP must itself have an 
energy density larger than the present DM density. 

To see if such production is sufhcient to give a dominant DM component, we need to know the 
yield of NLSPs after they have frozen out of the thermal plasma. For the neutralino NLSP yield, 
relevant processes include pair-annihilation and co- annihilation with the charginos, next-to-lightest 




(3.10) 
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Figure 5. The same as Fig. 4 but with the PQ scale /a = 5 x 10-* GeV. 



neutralinos and sleptons. For the stau NLSP, the yield is determined by the stau-stau annihilation 
and stau-neutralino co-annihilation processes. A typical relic abundance is 

r, = (i~io)xio-(^). (3^11, 

for a bino-dominated neutralino, and 

for the stau. Note that in the latter case the NTP can produce sufficient axino abundance to explain 
the whole DM density only for stau masses above 1.9 TeV, which may thermalize only at accordingly 
higher temperatures. 

These two choices for the NLSP were considered in the Constrained Minimal Supersymmetric 
Standard Model (CMSSM) in Ref. [4], for fa < 10^^ Qgy, for which even the stau lifetime is of order 
Is, or less. Recently, the case of stau NLSP, including four-body hadronic decays, was discussed in 
Ref. [8] also for larger values of fa- 
in conclusion, for neutralino NLSP, which decays mainly into an axino and a photon or a Z-boson, 
the NTP production is usually more efficient. For the stau NLSP, which can decay to an axino and a 
tau-lepton through a coupling of the type given in Eq. (3.6), the contribution is smaller, but can still 
be substantial. 
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Figure 6. The same as Fig. 4 but with DFSZ model used in Fig. 3. 



Regarding other NLSPs, colored relics (or even a wino-like neutralino if it is lighter than 1.8 TeV) 
usually remain in thermal equilibrium so long that their number density after freeze-out becomes 
negligible and therefore cannot produce any substantial axino population after freeze-out [37, 38]. 

4 Cosmological constraints 
4.1 The relic density of dark matter 

For the total axino DM relic density, we apply the 3cr range derived from WMAP-7 data [39] 

0.109 < r^a/i^ < 0.113. (4.1) 

This produces a stripe in the parameter space and also plays a role of the upper bound on the relic 
density when there are additional DM components, e.g. the axion. 

This can be seen in Figs. 4 and 5, where we present the reheating temperature versus the axino 
mass plane for /„ = 10^^ GeV and 5 x 10^ GeV, respectively. We have included both the thermal and 
the non-thermal production contributions of axinos, the latter assuming Ynlsp = (black solid) 10"^'' 
(green sohd) and 10~* (red dash), denoted also as (I), (II), and (III), respectively. A typical stau and 
neutralino yield after freezout will lie between (I) and (III). A correct relic density of axinos, in the 
range given by Eq. (4.1), corresponds to the thin bands between like curves. The parameter space 
above the curves is excluded as giving too much relic abundance. Similar figures in the DFSZ model 
are presented in Figs. 6 and 7, which can be compared to the figures in Ref. [47]. 
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Figure 7. The same as Fig. 6 but with the PQ scale /a = 5 x IQ-* GeV. 



4.2 Nucleosynthesis 

An injection of high energy electromagnetic and hadronic particles during or after Big Bang Nucle- 
osynthesis (BBN) epoch may disrupt the abundances of light elements. For axino DM, the lifetime of 
the NLSP, such as the neutralino or the stau, is typically around 1 sec, or less, and therefore constraints 
from the BBN are weak. However, for longer lifetimes such constraints become important [40, 41]. This 
leads to an upper bound on the decay products of the NLSP for a given NLSP lifetime. For the stau 
NLSP, a bound state with ^He severely constrains its lifetime to be less than roughly 5 x 10"^ sec [42] 
(although in specific cases with gravitino as DM, this can be up to an order of magnitude larger [43]). 
However for the parameters considered in our study, i.e. fa < 10^^ GeV, the BBN constraint can be 
avoided due to the small lifetime of the NLSP. For larger values of fa, on the other hand, non-trivial 
constraints arise, especially for the stau NLSP, as recently discussed in Ref. [8]. 

4.3 Structure formation 

The density perturbation due to axino population is suppressed at scales below their free-streaming 
length. When the axino mass is larger than ©(keV), thermally produced axinos become cold [2]. 
However, the non-thermal population of axinos from NLSP decays can still have a large velocity 
dispersion and can be too warm. Lyman-a [44] and reionization data [45] give a bound on the 
velocity of the WDM component and its fraction in the DM density. A recent analysis using the 
SDSS Lyman-a data [44] leads to an upper limit on the average velocity, (i')wdm < 0.013 km/s for 
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pure WDM, or otherwise in the case of mixed cold/warm DM the WDM fraction is constrained to be 
^wbm/^bm < 0.35 in the larger velocity region. 

The present velocity of thermally produced axinos is given by [46] 

vo 0.065 km/s ( . (4.2) 

Therefore the above Lyman-a data implies ma > 5 keV for TP axinos. This lower bound is marked 
in Figs. 4 and 5 with a vertical blue dashed line and an arrow. 

The free-streaming velocity of axinos produced non-thermally can be obtained from the lifetime 
of the NLSP and the mass relations [2, 4, 45]. For the bino-like neutralino NLSP with CaYY = 8/3, 
we find 

^ 0.4km/sf ( "^^ ) ( ] , (4.3) 

" ' VlMeV/ VlOOGeV/ V10"GeVy' ^ ' 



and for the stau NLSP, 
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For NTP axinos and for the neutralino NLSP we find a lower limit on the axino mass ma > 30 MeV 
and 1.5 MeV when fa = 10^^ GeV and 5 x 10^ GeV, respectively. These limits are marked in Figs. 4 and 
5 with a vertical blue dashed line and an arrow. For the stau NLSP the analogous lower bounds are 

> 150 MeV and 7.5 MeV, respectively. We stress that these bounds apply solely if the population 
of axinos produced through NTP is substantial (> 20 — 30%). 

In Figs. 8 and 9, we show contours of the reheating temperature in the plane spanned by 
the NLSP and the axino mass. ^° In Fig. 8 we have fixed fa = 10^^ GeV and assumed Ynlsp — 
10~^^ (?7iNLSp/100 GeV), which is a typical value for bino-like neutralino NLSP. The cyan wedge in 
the upper right-hand corner is excluded by the overdensity of DM, while in the red wedge below it the 
axino is not the LSP. In Fig. 9 instead Inlsp — 10"^'' (mNLSp/100 GeV) has been assumed, typical of 
stau NLSP. So long as the TP dominates, the curves of constant Tr remain vertical and practically 
independent of ma, while as soon as the NTP becomes important, ma dependence arises leading to 
non- vertical curves. For example, with a bino-like neutralino NLSP of 100 GeV, as in Fig. 8, the NTP 
contributes only up to 10 % of the axino LSP DM density. For stau NLSP with small abundance 
Y < 10"^^, as given in Eq. (3.12), NTP is always subdominant so that the contours remain vertical. 
In this case the bounds coming from the free-streaming velocity of the NTP axinos are absent. 



5 Conclusion 

We have performed an updated analysis of the relic axino production, taking into account some new 
calculations that have appeared after the initial study [1, 2], especially [5, 6] for the thermal part, and 
compared them with our own results to explore the question of uncertainty and model-dependence. 

We have found that the uncertainty has not really decreased after the latest calculation [6]. This 
is probably not surprising since the QCD coupling is large and the convergence of the perturbative 

^"Similar figures are shown in Ref. [8]. 



- 22 - 



Choi, Covi, Kim, Roszkowski '11 

I I I I III il\ I I I I lllj I I I Mil 




Figure 8. Contours of constant reheating temperature in the NLSP-axino mass plane. Here we have assumed 
"Knlsp = 10"^^ (rriNLSp/lOOGeV), typical of neutralino NLSP, and taken fa = lO" GeV. The cyan wedge in 
the upper right-hand corner is excluded by the overdensity of DM, while in the red wedge below it the axino 
is not the LSP. 

series is quite slow. Comparing the different results, we estimate the uncertainty in the relic density 
of axinos produced thermally to be still of order a factor of 10 or so at Tr ~ 10'' GeV, and to that one 
has to add also some possible (unknown) contributions due to non-perturbative effects. Our result lies 
above both estimates given in Refs. [5, 6], and this seems natural since we included subleading terms 
in m^g/s, which do indeed increase the cross-sections for single channels ensuring their positivity in 
the whole range of integration, even in the limit of very large gauge coupling. 

Regarding the model dependence, our conclusions are more optimistic: in the KSVZ-type the QCD 
anomaly term strongly dominates the axino thermal production mechanism, apart from the case of 
small reheating temperatures where sparticle decay contributions start playing the dominant role. The 
inclusion of the additional anomaly couplings is completely negligible, apart from unnaturally large 
values of the coefficients Caww and CaVY- Instead for DFSZ-typc models, the Yukawa interaction 
dominates around the weak scale and can give the main production mechanism of axinos making it 
independent of the reheating temperature. Therefore the axino abundance is free from the uncertainty 
in the method used for the IR-divergence. At large temperatures, it is the EW anomaly term that 
dominates, giving a lower abundance than in the KSVZ case. 

For both models, also the non-thermal production via NLSP decay can produce the required axino 
DM density, if the NLSP decouples with a sufficiently large abundance. But for this mechanism to 
dominate, the reheating temperature has to be very low and the axino and NLSP masses not too 



- 23 - 



Choi, Covi, Kim, Roszkowski '11 




10' 

■axmo(GeV) 



Figure 9. The same as Fig. 8 except for Ynlsp = 10 (ttinlsp/IOO GeV), typical of stau NLSP. 



hierarchical. We find that, interestingly enough, a light bino NLSP decaying out of equilibrium can 
still produce the whole DM density at the cost of a very low reheating temperature (but sufficiently 
high for NLSP thermalization). For the stau case instead, it is quite unlikely that the NTP can 
dominate, unless the stau NLSP yield after freeze-out is unusually large. 

During the final stages of completion of this work, the analysis Ref. [47] appeared which discusses 
in details axino couplings and finds a non-trivial momentum dependence in the one-particle irreducible 
one-loop axino-gluon-gluino couplings. The coefficient Cip/ of these interactions is suppressed when 
the external particle momentum is much larger than the mass M of the PQ-charged fermions in the 
loop. Due to this effect, the authors claim a suppression of order Af^/T^ of the axino production 
from the dimension-5 operators for the DFSZ case and for extremely small KSVZ quark masses (with 
Yukawa couplings less than 10~^, i.e. for the heavy quark mass less than 10^ GeV for fa — 10^^ GeV). 

We investigated such suppression by inserting their Cip/ coupling in the relevant diagrams and 
we obtain different suppressions depending on the type of Feynman graph and a strong dependence 
on the IR regulator contained there, in most cases the gluon thermal mass. In particular, we find 
no suppression at all for the t-channel gluon exchange for vanishing gluon mass. Since graphs with 
the one loop Cip/ coupling and a gluon thermal mass insertion arise at lowest order at two loops, 
probably a full investigation of the two-loop diagrams in thermal field theory is needed to resolve this 
issue. 

Note in any case that even without suppression, the DFSZ axino production is dominated by the 
decay term up to temperatures of the order of 10^~^ GeV. For the KVSZ case, we show in Fig. 10 as 
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Figure 10. The same as Fig. 4 but showing as well in magenta the yield suppression found in Rcf. [47] for 
different values of the heavy quark masses. 

violet lines how the yield changes according to Ref. [47] for small heavy quark masses, mg — 10^ GeV 
and mq — 10^ GeV. The Cold DM axino, on vifhich our present study is based, is practically not 
affected. For large fermion masses in the loop, M > T, our anomalous couplings coincide fully with 
the Cipj in Ref. [47] and our results are in perfect agreement. 
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A Appendix 



In this Appendix, we present the computation of the axino rehc density in detail. The time evolution 
of the axino number density, , is described by the Boltzmann equation: 

+ SHria = E(o-(« +j^a-\ )vrci)ninj + ^(r(i ^ 5 H ))ni. (^.l) 



Here H is the Hubble parameter, H{T) ~ \J {ir'^g*)/ (90Mp)T^, where is the number of effective 
massless degrees of freedom. The first term in the r.h.s. of Eq. (A.l) is the axino production from the 
scattering process of particles i and j in the thermal bath with cross section cr(i + j — > 5 + • • • ) and 
relative velocity Wrei- ni is the «-th particle's number density in the thermal bath. The second term 
in the r.h.s. of Eq. (A.l) is the axino production from the decay of i particle which are in thermal 
bath with the decay rate r(i — )• a + • • • ). (• • •) denotes the thermal averaging including averaging over 
initial spins and summing over the final spins. Here we neglect the inverse processes since the axinos 
are decoupled from the thermal bath and their number density is very small. 
Using the axino yield defined as 

F^^, (A.2) 

where s = {2Tr'^ / 45) gs*T^ is the entropy density, and normally gs* = g* in the early Universe, the 
solution of the Boltzmann equation is 



(A.3) 



where 



Jo sHT 



(A.4) 



sHT 



The explicit formulae for the thermal average are given in Ref. [48]. 

The relevant 2-body scattering cross sections are summarized in table 2, keeping in mind that the 
physical cross-sections are Eq. (3.2) 

Then tT„(s) is given from the matrix element, M„, by 

^n{s) = ^4 f Wn?dt. (A.6) 



Among the processes, B,F,G and H have an infrared (IR) divergence due to the massless gluon exchange 
in the t- or M-channel. To regularize this IR divergence, we introduce, only in the t or u-channel gluon 
propagators, an effective thermal gluon mass m^fj = g^T^, which is generated by plasma effects [2, 49]. 
Then we obtain finite and always positive cross sections. For example, the squared matrix element for 
process B with using massless gluon propagator is 

\Ms?^-'^^'"\^''^'"\r^^\\ (A.7) 
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However introducing effective thermal gluon mass in the t- or w-channcl, it changes to 



which is positive definite by definition. Then crs(s) is obtained after integration over t as given in the 
Eq. (A. 6). Here we list the full cross-sections (t„ for these processes. 



1. Process B 

abc\2 



jfz' ^ ^ b (3"^off - 5sm2ff - 7s^) + (-3™^^ + dsmt^ + I2s^ml^ + 45^) log ((m^g + s)/™^ )] 

(A.9) 

2. Process F: We include the factor 1/2 for the identical initial particles, 

12s2 (27712^. + s) L V cff off ; ^^^^Q^ 

+12 {2ml^ + Qsmls + As^mls + s^) log {{ml^ + s)/mls)] . 



3. Process G: 



4. Process H: 



^ [-2. + {2ml^ + s) log {{ml^ + s)/ml,)] . (A.ll) 



-s (3m4g + 9sm2g + 7s^) + (3m^ff + ILsm*^ + 125^^2^ + 45^) log ((m^j^ + s)/mlg)] 



16s2 {mlg + s) 

(A.12) 

These formulas give back the expressions in Table 2, in the limit of large s/rrics. We mention 
here that those expressions in Table 2 provide unphysical result for small s/rrics and give negative 
contribution to the thermal averages. In fact substituting s ^T^, m^g ~ g^T^, one finds e.g. for the 
B process 

as cx log[l/g2(r)] - I < for g^T) > 0.17 (A.13) 

and this holds for all the interesting region up to the GUT scale where g'^ ^ 0.5. So even if the total 
thermal cross-section may result positive, it is lower than the real one due to the negative contribution 
of these IR-divergent channels. On the other hand the full expressions above from Eq. (A.9) to 
Eq. (A.12) are positive for any value of s. We can see that easily in the limit of large mlg/s because 
those reduce to: 

1. Process B: 

3|jQbc|2 

32 

2. Process F: 



(A.14) 



' '0{^] (A.15) 



24 



'eft 
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3. Process G: 




+ 



( 



) 



(A.16) 



24 mt, 



eff 



6 

'eff 



4. Process H: 



32 



+ 



TO- 



S 



2 

'off 



) 



(A.17) 



so they have the correct asymptotic behaviour corresponding to screening and physical decouphng of 
the intermediate gluon channels for large gluon mass. So with these formulas the single cross-sections 
are always positive and the thermal average larger than in the previous estimates. 

In our numerical computations we use the full formulae above for the IR divergent processes, while 
for the other processes, A, C, D, E, I, J, we use the expressions are given in Table 2. 
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